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ABSTRACT 

In this letter we report the possible existence of a second planet in the tran- 
siting extrasolar planet system HD 17156 and its interactive dynamics with the 
previously known planet. The analysis is achieved through the pofp optimiza- 
tion software which is based on a full integration of the system's multiple-body 
Newtonian equations of motion. The two-planet solution yields a significantly 
improved fit to the previously published radial velocities. The two planets are 
strongly interacting and exchange angular momentum in a 5:1 mean motion res- 
onance, yet remain stable as they mutually excite orbital eccentricities and peri- 
astron advances. 

Subject headings: planetary systems - celestial mechanics - gravitation - insta- 
bilities - methods: N-body simulations - methods: numerical 



1. Introduction 



By virtue of its unusual characteristics, HD 17156b is one of the most valuable extrasolar 
planets for understanding planet formation and orbital dynamics. Discovered via the Doppler 
technique by the N2K consortium (Fischer et al. 2007), the planet was found to transit its 
host star by the TransitSearch.org collaboration (Barbieri et al. 2007). Additional 
transit observations and refined system parameters are given by Gillon et al. (2008), Narita 
et al. (2008), and Irwin et al. (2008). The planet radius and mass are approximately 1 Rj up 
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and 3 Mj up , giving a density of ~ 3.5-4.0 g cm -3 , nearly three times that of Jupiter. More 
interesting however are the planet's orbital characteristics: (i) The 21.2 d orbital period is 
roughly a factor 7 times longer than mostfj] other transiting planets; (ii) While most hot 
Jupiter planets have low eccentricities (the majority consistent with zero), HD 17156b has 
a very high eccentricity of e=0.67. Although the high eccentricity coupled with the small 
semimajor axis (0.16 AU) does not necessarily require the presence of a third body perturbing 
the planet's orbit (Gillon et al. 2008), it suggests that such a body may be present. In this 
work we examine the published radial velocities of Fischer et al. (2007) and the transit times 
as given in Irwin et al. (2008) using a Newtonian (not Keplerian) 3-body integrator and 
conclude that a second planet is not only possible, but probable. Notably, while the two 
planets exhibit strong dynamical interaction leading to large eccentricities and periastron 
advances, they remain stable via an elegant exchange of orbital angular momentum. 



2. Initial Analysis of Radial Velocities 

Our initial investigation of the HD 17156 system began with a 1-planet fit to the Keck 
and Subaru radial velocity data from Fischer et al. (2007). We omitted the radial velocity 
data of Narita et al. (2008) since these data are influenced by the Rossiter-McLaughlin effect, 
but their inclusion has essentially no effect on the results. We reproduced the parameters of 
the planet HD 17156b as reported in Fischer et al. (2007), though we found that an additional 
0.14 ms _1 offset added to the Subaru data set gave a slightly better fit. We computed a 
power spectrum of the residuals of the radial velocities after removing the one-planet fit. A 
Monte Carlo technique was used to allow us to assess the significance of any peaks. The 
power spectrum indicated a possible peak with a period of approximately 115 days. This 
prompted us to undertake a much more thorough and realistic 3-body investigation described 
below. In addition to the power spectra, we computed Keplerian periodograms: for a given 
period the parameters of a Keplerian fit to the radial velocities were optimized (K, e, u>, T , 
7), i.e. the Lehmann-Filhes equation (see Hilditch 2001 for example). While a traditional 
Fourier power spectrum has no free parameters, it uses sines and cosines as basis functions. 
However, an eccentric orbit is far from being sinusoidal and as a consequence a Fourier power 
spectrum will require power at many periods to match the radial velocities. In contrast, by 
using the Lehmann-Filhes equation as the basis function we obtain an optimal periodogram. 
Furthermore, unlike a power spectrum, fitting a Keplerian orbit at each trial period allows 
one to use the error bars on each observation as weights. Thus no Monte Carlo or bootstrap 
sampling is necessary to assess the quality of the fit: the periodogram provides a \ 2 versus 
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period directly. Such periodograms clearly showed the 21.2 d signal from planet b, and the 
residuals again indicated the presence of a possible second planet with a period of ~106-116 d. 

3. HD 17156c 

To investigate characteristics of a possible second planet we used the Planetary Orbit 
Fitting Process (pofp), an optimization software written and operated in matlab — see 
Windmiller, Short & Orosz (2007) and Short, Windmiller & Orosz (2008). The pofp provides 
multi-body solutions based on a full integration of the Newtonian equations of motion. 

In addition to the radial velocity data, we used pofp to fit the 4 transit times as 
listed in Irwin et al. (2008). The inclusion of the transit times was crucial because it not 
only greatly improved the period of the known planet b, but the precise timings strongly 
constrain the orbital configuration and interaction between the planets. Since the system 
transits, we assumed the second planet's orbit is edge-on and co-planar with the first. A 
careful reading of Fischer et al. (2007) led us to use caution with the Subaru radial velocities, 
thus we considered two data set combinations: The first consisted of the 24 high-quality radial 
velocities from Keck plus the 4 observed transit times, and the second consisted of the 33 
combined Keck and Subaru radial velocities plus the 4 transit times. Following Fischer et 
al. (2007), we set the stellar mass to 1.2 M Q and the stellar jitter to 3 ms" 1 . Using pofp wc 
fit the two data sets with a single planet model giving two solutions denoted here by POFPl 
K+TT (Keck plus transit times) and POFPl K+S+TT (Keck, Subaru, plus transit times). 
This provided a baseline set of models to compare against when considering a two-planet 
solution. The orbital parameters are given in Table 1, along with the Fischer et al. (2007) 
and Irwin et al. (2008) solutions. 

We then fit using the Keck and Subaru radial velocities plus transit times with a two- 
planet model, giving the solution POFP2. Because the POFP2 solution is fully Newtonian 
and the planets interact, the tabulated parameters are valid only at a specific time, which 
is taken to be the time of the first Keck radial velocity measurement. The fitness for the 
1-planet and 2-planet solutions is given in Table 2. Note that in this table only the following 
are optimized: the 1-planet model using Keck radial velocities, the 1-planet model using 
Keck+ Subaru velocities, and the 2-planet model using Keck+Subaru velocities; in all cases 
the transit times were used. These optimized values are shown in bold-face; the other 
entries result when these models are evaluated and matched against the listed data set and 
are given for comparison. Some reduced x 2 values are less than 1.0, suggesting that the 
3 ms -1 jitter estimate of Fischer et al. (2007) is slightly too large. The key point of Table 2 
is the following: including two planets significantly improves the fit in all cases (Keck only, 
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Keck+TT, Keck+Subaru, Keck+Subaru+TT). 

In our 2-planet solution we found that the inner planet maintains an almost constant 
semi-major axis of 0.1596 AU, a slowly oscillating eccentricity having a period of 141.1 yr 
with e between 0.665 and 0.670, and a full rotation of its line of apsides every 33,050 yr. 
Thus, there are two long term cycles of approximately 141.1 yr and 33,050 yr. Meanwhile, the 
outer planet HD 17156c exhibits an oscillation in eccentricity with a pronounced amplitude, 
varying between e = 0.10 and e = 0.50, in a near-mirror image of the eccentricity oscillation 
of planet b (see Fig. 1). These complementary oscillations of eccentricities provide a very 
interesting example of planetary angular momentum interchange. The semi-major axis of 
planet c varies between ~ 0.46 and 0.49 AU and the position of its periastron completes a 
141.1 yr cycle in step with the change in the eccentricities. Furthermore, for each 141.1 yr 
cycle, there will be an orbit of maximum eccentricity for planet c. Over 33,050 yr the 
periastron of this maximum eccentricity orbit will advance slowly and synchronously with 
the precession of planet b's periastron. As a result of this resonance coupling between the 
orbits, the stability of this system is maintained by keeping an approximately 0.14 AU buffer 
between the two planets' orbital paths. This is illustrated in Fig. 2, showing the change in 
the orbital configuration caused by planetary interaction. Over a span of half the 141.1 yr 
cycle we see planet c's eccentricity change dramatically. The lower left panel shows a much 
larger advance in time and is a snapshot of a maximum eccentric orbit of planet c; notice 
the simultaneous precession of both planets' periastra at the longer 33,050 yr cycle. 

In addition to any goodness of fit analysis, the two-planet model must show long-term 
stability. We have used the symplectic integration packaged HNBody of Rauch and Hamil- 
ton (2002) to do a long term integration of this system. Over a time interval of 100 million 
days (~274,000 yr), we have verified stability and that the planets exhibit a 5:1 mean motion 
resonance. Following the lead of Barnes (2006), we plotted et,e c sin(Ac<j) vs. e;,e c cos(Acj) in 
Fig. 3, where Au is the difference in the longitude of the planets' periastra. In this figure, 
if the angle between the periastra was constant, all the points would be co-linear, extending 
radially outward from the origin. The closed loop shape arises because the angle between 
periastra grows and spans a full 2n every 141.1 yrs; there is no libration. From inspection 
of the figure, the maximum of ef,e c occurs when Au = 0, and since e& ~ 0.67, the orbits of 
planet c reach their maximum eccentricity when the periastra are aligned. Hence the specific 
cases illustrated in Fig. 2 are in fact typical. 

The l-o" confidence limits for the parameters of the 2-planet solution were calculated 
following Press et al. (1986) and are given in Table 3. Note that we quote the parameter 



2 http:/ /janus. astro. umd.edu/HNBody/ 



- 5 - 



values to more significant digits than is warranted by the uncertainties so that others can 
exactly reproduce our 2-planet pofp2 numerical solution. Also note that while the short- 
term orbital characteristics of planet c can be determined, the timescales of the much longer 
precession cycles are poorly constrained by the data (e.g. an earlier POFP2 solution yielded 
long term cycles of 212.7 yr and 37,700 yr). 

The question remains, does HD 17156c exist? Noting the small number of observations, 
the short time span of these observations, and the small mass of the inferred second planet, 
caution must be taken. Using the F statistic on the fitness values from Table 2 to test if 
the extra parameters of the two-planet model are warranted, we find that the probability for 
including those additional parameters is greater than 99.9%. While the mass of HD 17156c 
is poorly determined, it is not consistent with zero, further supporting the F statistic result. 
The two-planet model also is physically realizable as a stable system. These facts led to our 
characterization of HD 17156c as not only possible, but probable. Obviously, to resolve the 
remaining uncertainty, additional observations are required. In §4, we point out observations 
that would allow distinguishing between the one and two-planet models. 



4. Discussion 

Additional radial velocity and transit time data would greatly help confirm and con- 
strain specific parameters of the second planet. However, the general dynamics of the 3-body 
interactions are already well determined. Using the current observations and our two-planet 
pofp model we can make several observable predictions. The divergence between the ve- 
locities derived from the single-planet reproduction and the two-planet solutions is shown in 
Fig. 4. In about 2 years, the differences in radial velocity predictions between the two models 
become larger than 10 ms" 1 at the peaks of the graphs, occurring at the time of periastron 
of planet b. In particular, for the epoch of planet b's periastron at HJD 2455266.04 (UT 
2010 Mar 10), a difference of 18 ms' 1 is predicted. 

Another observable is planet b's time of transit. While this is slowly advancing on the 
33,050 yr cycle, there is a much larger and rapid modulation in planet b's instantaneous 
period, i.e. the time between sequential transits. The near-future predicted instantaneous 
period is shown in Fig. 5. Planet-planet interactions, especially when near apastron of planet 
b, result in sharp changes in the instantaneous period. The average time interval between 
successive conjunctions is ~26.2 d which is ~24% longer than planet b's period. Thus 
conjunctions close to the apastron of planet b occur roughly every ~5 orbits. The exchange 
of angular momentum at these times causes the jumps in period of planet b shown in Fig. 5. 
A longer ~1000 d cycle results from the longitude of planet c shifting by -6.75 degrees during 
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each successive apastral conjunctions. A predicted O-C diagram using the mean period at the 
current epoch is shown in the lower panel of Fig. 5. Unlike the instantaneous period changes, 
an O-C diagram includes the sum of the effects of the period changes since the fiducial epoch, 
and thus is akin to the integral of the instantaneous period. The deviations from a linear 
ephemeris are quite large and thus monitoring planet b's transit times should reveal the 
presence of planet c and be very helpful in further refinement of the system's parameters. 
Finally, planet c may in principle also be detected via the transit method. Unfortunately, 
however, planet c is not expected to transit the host star if its orbit is coplanar with that of 
planet b (e.g. i = 86.5°, Irwin et al. 2008). 

Whereas the parameters of the probable planet c are poorly constrained, the qualitative 
dynamical aspects of such a system are of great interest. This case provides an example of a 
strong dynamical interaction between a high eccentricity "hot Jupiter" planet and a second 
planet found further away from the star. The predicted near-term observables should provide 
guidance to future observations of this system and may provide a possible explanation to 
any variations seen in observed transit times. 

We thank the anonymous referee for suggestions that helped improve this work. We 
thank Rory Barnes of the Lunar and Planetary Laboratory of the University of Arizona for 
his timely and thoughtful comments on the dynamics of the two-planet model. 
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Table 1. HD 17156 Planet Parameters 



Solution 


Epoch a 


P 


T b 


e 


LU 


K 


a 


Mass d 




(HJD-2450000) 


(days) 


(HJD-2450000) 




(deg) 


(ms _1 ) 


(AU) 


(M Jup ) 


Fischer et al. 2007 




21.2 


3738.529 


0.67 


121 


275 


0.15 


3.12 


Irwin ct al. 2007 (K+S+TT) 




21.2169 


3738.605 


0.670 


121.3 


273.8 


0.160 


3.13 


pofpI K+TT 




21.2167 


3738.614 


0.670 


121.4 


273.8 


0.160 


3.12 


pofpI K+S+TT 




21.2168 


3738.614 


0.670 


121.6 


273.3 


0.160 


3.13 


POFP2 Planet b 


3746.7582 


21.2144 


3738.593 


0.670 


120.9 


273.0 


0.160 


3.13 


and Planet c 


3746.7582 


111.394 


3736.880 


0.136 


273.6 


2.6 


0.481 


0.07 



a Time at which tabulated parameters are valid for the POFP2 solutions 
b Time of periastron passage 

c Assuming a stellar mass of 1.2 Mq from Fischer et al. 2007 



Table 2. Fitness a of Solutions 13 





Number of 


rms 


Keck RV only 


Keck+TT 


Kcck+Subaru RV only 


Kcck+Subaru+TT 


Solution 


parameters 


(ms _1 ) 


N=24 


N=28 


N=33 


N=37 


Fischer et. al 2007 




3.97 






(1.37; 1.08) c 




POFPl K+TT 


5 


3.65 


26.47 (1.39) 


27.45 (1.19) 


30.41 (1.13) 


33.12 (1.07) 


pofpI K+S+TT 


6 


3.75 


26.49 (1.39) 


29.99 (1.30) 


30.48 (1.13) 


31.51 (1.02) 


POFP2 K+S+TT 


11 


3.14 


16.56 (1.18) 


17.90 (0.99) 


19.96 (0.91) 


21.38 (0.82) 



a Reduced x 2 values are in parentheses. 

b Optimized solutions are shown in bold; other values show the solution evaluated against the specific data set. 
c Note Fischer et al. (2007) give two values for the reduced \ 2 ■ 



Table 3. HD 17156 POFP2 Error Estimates 



Planet 


P 


To 


e 


LU 


K 


a 


Mass 




(days) 


(HJD) 




(deg) 


(ms- 1 ) 


(AU) 


(M Jup ) 


+ 1(7 


+0.0031 


+0.05 


+0.003 


+0.14 


+1.77 


+0.0001 


+0.02 


Planet b 


21.2144 


3738.593 


0.669614 


120.89 


272.987 


0.159505 


3.125 


-la 


-0.0025 


-0.05 


-0.0005 


-0.03 


-1.75 


-0.0001 


-0.02 


+ 1(7 


+0.43 


+0.58 


+0.0029 


+0.17 


+0.79 


+0.001 


+0.021 


Planet c 


111.3937 


3736.880 


0.136492 


273.57 


2.580 


0.481478 


0.068 


-1(7 


-0.09 


-1.61 


-0.0012 


-0.17 


-0.71 


-0.0003 


-0.019 



a Full precision is given so the initial conditions of the solution are exactly specified. 
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Fig. 1. — The cyclical nature of the planets' orbital properties is seen in this POFP2 inte- 
gration spanning 10 6 d. The panels show, from top to bottom, the eccentricity of planet b 
(eft), the eccentricity of planet c (e c ), and the longitude of periastron of planet c (u> c ). The 
141.1 yr (51,550 d) periodicity is manifest in the small modulation of eft (upper panel) and 
the much larger changes in e c and u c . The anti-phasing of the eccentricities is a result of the 
exchange of angular momentum between the planets. 
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Fig. 2.— The orbital configuration of the POFP2 two-planet solution is shown at four 
separate epochs. Clockwise from the upper left, the panels show the advance of planet c's 
cycle (141.1 yr) at approximately 1/4 cycle intervals. The lower left panel illustrates the 
precession of planet b's orbit, on its much longer cycle of 33,050 yr. The cross marks the 
periastron of planet c's orbit, and the arrow in the lower right corner marks the direction to 
the observer. 
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Fig. 3. — The resonant coupling of the eccentricities of the two planets is shown, following 
the suggestion of Barnes (2006) of plotting the product of the instantaneous eccentricities 
times the sine and cosine of the angle between the two planets' periastra (Au). The orbits 
were sampled every 50,000 days over the 10 8 d interval of the HNBody integration. Since 
the eccentricity of planet b is almost constant, the figure shows that the orbit of maximum 
eccentricity of planet c occurs when the periastra of the two planets are approximately 
aligned, i.e., along the abscissa where Au ~ 0. 
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Fig. 4. — The difference between the radial velocity predictions of the pofp one-planet and 
two-planet solutions is shown. The bottom panel is an enlargement of a 50-day interval from 
the upper panel. The largest differences (peaks in the lower panel) occur around the time of 
planet b's periastron, and become readily measurable around HJD 2455200 (~ 2010 Jan). 
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Fig. 5.— Upper panel: The instantaneous period of planet b in the POFP2 two-planet 
solution plotted against time. Angular momentum exchange between the planets results 
in jumps in the orbital period. Changes of up to ~3 min in the transit-to-transit period 
are predicted. Lower panel: A predicted O-C diagram for the times of transit of planet b, 
with the transits times from Irwin et al. (2007) superposed, using the ephemeris T(E) = 
2453738.32783 + 21.216159 E (HJD). 



